
; This NCL script takes an existing SOM forcing file and adds heat sources and sinks 
; the anomaly is described as either annual average, which is put into every month of the qdp variable, or only as an anomaly in some month
; used normal pop_frc.gx1v6.100513.nc file (used by Kate; downloaded from https://svn-ccsm-inputdata.cgd.ucar.edu/trunk/inputdata/ocn/docn7/SOM/) 
; as the base state

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 

case = "annual avg qdp"

begin

 old_input = "mk_qflux_forcing_files/b.e104.B_1850_CN.f19_g16.qflx_clim.nc"   ; selfmade control file
 f1 = addfile(old_input,"r")

 area  = f1->area
 maskr = f1->mask
 xc    = f1->xc
 yc    = f1->yc
 time  = f1->time
 S     = f1->S
 T     = f1->T
 U     = f1->U
 V     = f1->V
 dhdx  = f1->dhdx
 dhdy  = f1->dhdy

 qdp0  = f1->qdp
 hblt  = f1->hblt 
;--------------------------------------------------------------------- for weighting
;dir_path_in   = "/net/meso/climphys4/mariaru/cesm104/abrupt4x/prepared/"
;latlonO       = addfile(dir_path_in+"grid_kmt.nc","r")
dir_path_in   = "/home/divanova/regridding/grids/gx1v6/"
;dir_path_in   = "/home/divanova/Grids/gx1v6/"
latlonO       = addfile(dir_path_in+"grid_gx1v6.nc","r")

TLAT          = latlonO->TLAT
TLONG         = latlonO->TLONG
TAREA         = doubletofloat(latlonO->TAREA)/10000 ; from cm^2 to m^2
swgt          = sum(TAREA)
print("swgt is "+swgt)
delete(latlonO)

qdp0!0 = "month"
qdp0!1 = "lat"
qdp0!2 = "lon"
qdp0@lon2d = TLONG
qdp0@lat2d = TLAT
copy_VarMeta(qdp0,hblt)
;printVarSummary(TLAT)
;print(TLAT(:,285))
;print(TLAT(:,35))


;-------------------------------------------------------------------- put in Q fluxes 
qdp_anom = qdp0 ;init
qdp_anom = 0

pi  = 3.1415
; test combination Ind_med: neg Indian mag = 23 (as t3 before) 
; test combination Ind_stg: neg Indian mag = 47 (as t4 before) 
; test combination Ind_low: neg Indian mag = 23/2 (half t3 before) 

; magnitudes from tests before
; test combination t1: pos. Indian mag = 47, neg. subtrop mag = 61.05
; test combination t2: pos. Indian mag = 23, neg. subtrop mag = 30.4
; test combination t3: neg. Indian mag = 23, pos. subtrop mag = 31
; test combination t4: neg. Indian mag = 47, pos. subtrop mag = 61.7

;------------------ tropical Indian Ocean
mag = 47
do i=70,130-1
low = ind_nearest_coord(-20,TLAT(:,i),0)   ; Eq at that point on the kmt grid 
up  = ind_nearest_coord(20,TLAT(:,i),0)  ; 10N at that point on the kmt grid 
gradient = fspan(-pi/2,pi/2,up-low+1)
qdp_anom(:,low:up,i)  = conform(qdp_anom(:,low:up,i),cos(gradient)*mag,(/1/))
;qdp_anom(:,low:up,i)  = conform(qdp_anom(:,low:up,i),-cos(gradient)*mag,(/1/))
delete(gradient)
end do

;------------------------------- prep writing out
qdp = qdp0 + qdp_anom
copy_VarMeta(qdp0,qdp)

 ;test and write out global q flux
 Q = month_to_annual(qdp,1)
 Qsum = sum(where(Q(0,:,:).eq.-999,9.96921e+36,Q(0,:,:)*TAREA))/swgt
 print("total q flux is = "+Qsum)

 ;------------------------------ write out
; system("rm -f b.e104.B_1850_CN.f19_g16.qflx_clim_ind_low.nc")
; fout = addfile("b.e104.B_1850_CN.f19_g16.qflx_clim_ind_low.nc","c")
; setfileoption(fout,"DefineMode",True)

; fileAtt = True
; fileAtt@title = "annual mean ocean forcing from POP output"
; fileAtt@conventions = "CCSM data model domain description"
; fileAtt@source = "pop_frc.ncl"
; fileAtt@description = "Input data for DOCN7 mixed layer model from " + case
; fileAtt@note1 = "fields computed from disturbing the SOM forcing file b.e104.B_1850_CN.f19_g16.qflx_clim.nc"
; fileAtt@note2 = "all fields interpolated to T-grid"
; fileAtt@note3 = "qdp is computed from depth summed ocean column"
; fileAtt@author = "M. Rugenstein"
; fileAtt@calendar = "standard"
; fileAtt@comment = "This data is on the displaced pole grid gx1v6"
; fileAtt@creation_date = systemfunc("date")
; fileattdef(fout,fileAtt)
 
; fout->area = area
; fout->mask = maskr
; fout->xc = xc
; fout->yc = yc
; fout->time = time
; fout->S = S
; fout->T = T
; fout->U = U
; fout->V = V
; fout->dhdx = dhdx
; fout->dhdy = dhdy
; fout->hblt = hblt
; fout->qdp = qdp


qdp_ann  = qdp0(0,:,:) ;init
qdp_ann  = month_to_annual(qdp0,1)
qdp_ann1 = qdp0(0,:,:) ;init
qdp_ann1 = month_to_annual(qdp,1)
minus = qdp0(0,:,:);init
minus = qdp_ann1 - qdp_ann
;============================================ plotting
;       Sahel zone shape file
 dir="/home/divanova/SahelProject/topo/sahel_zone/"
 filename="sahelzone.shp"


wks_type = "png"
wks  = gsn_open_wks(wks_type,"qdp_t4.png")
  plres = True
  plres@gsLineColor = "red"
  plres@gsLineThicknessF =5

res                             = True
res@gsnDraw                     = False
res@gsnFrame                    = False
;res@vpWidthF                    = 0.8
;res@vpHeightF                   = 0.4
res@cnFillOn                    = True
res@cnInfoLabelOn               = False
res@cnLinesOn                   = False
res@cnLineLabelsOn              = False
res@gsnAddCyclic                = True
;res@gsnLeftStringFontHeightF    = 0.025
;res@gsnRightStringFontHeightF   = 0.025
;res@tmXBLabelFontHeightF        = 0.024
;res@tmYLLabelFontHeightF        = 0.024
;res@lbLabelBarOn                = False    ; no individual label bars
res@mpPerimOn   = True        ; Turn on map perimeter.
gsn_define_colormap(wks,"MPL_RdBu")
gsn_reverse_colormap(wks)
res@cnLevelSelectionMode        = "ManualLevels"
res@cnMinLevelValF              = -70
res@cnMaxLevelValF              = 70
res@cnLevelSpacingF             = 2.5
;res@tmXBLabelsOn                = False
;res@tmYLLabelsOn                = False
res@gsnLeftString               = "q flux anomaly, global = "+Qsum+" W/m~S~2~N~"
res@gsnRightString              = " "

res@gsnMajorLatSpacing = 30              ; change maj lat tm spacing
res@gsnMajorLonSpacing = 30              ; change maj lon tm spacing
res@tmXBMinorOn        = False           ; no lon minor tickmarks
res@tmYLMinorOn        = False           ; no lat minor tickmarks
res@tmYROn       = False           ; no lat minor tickmarks
res@tmXTOn       = False           ; no lat minor tickmarks

plot = gsn_csm_contour_map(wks,minus,res)
dum  = gsn_add_shapefile_polylines(wks,plot,dir+filename,plres)

draw(plot)
frame(wks)
;delete(res@cnLevelSpacingF)
res@gsnLeftString               = "annual average qdp without anomaly"
plot = gsn_csm_contour_map(wks,qdp_ann,res)
draw(plot)
frame(wks)
res@gsnLeftString               = "annual average qdp with anomaly"
plot = gsn_csm_contour_map(wks,qdp_ann1,res)
draw(plot)
frame(wks)

;delete([/res@cnLevelSelectionMode,res@cnMinLevelValF,res@cnMaxLevelValF/])
;res@gsnLeftString               = "hblt"
;plot = gsn_csm_contour_map(wks,hblt(0,:,:),res)
;draw(plot)
;frame(wks)

end
